library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
[30m── [1mAttaching packages[22m ────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──[39m
[30m[32m✓[30m [34mggplot2[30m 3.3.0 [32m✓[30m [34mpurrr [30m 0.3.3
[32m✓[30m [34mtibble [30m 3.0.0 [32m✓[30m [34mdplyr [30m 0.8.5
[32m✓[30m [34mtidyr [30m 1.0.2 [32m✓[30m [34mstringr[30m 1.4.0
[32m✓[30m [34mreadr [30m 1.3.1 [32m✓[30m [34mforcats[30m 0.5.0[39m
[30m── [1mConflicts[22m ───────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(cowplot)
********************************************************
Note: As of version 1.0.0, cowplot does not change the
default ggplot2 theme anymore. To recover the previous
behavior, execute:
theme_set(theme_cowplot())
********************************************************
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:cowplot’:
stamp
The following objects are masked from ‘package:dplyr’:
intersect, setdiff, union
The following objects are masked from ‘package:base’:
date, intersect, setdiff, union
library(mgcv)
Loading required package: nlme
Attaching package: ‘nlme’
The following object is masked from ‘package:dplyr’:
collapse
This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
bes<- read_csv("dataOut/binned_EachSize.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
bds <- read_csv("dataOut/binned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
ues <- read_csv("dataOut/unbinned_EachSize.csv")
Parsed with column specification:
cols(
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m,
depth = [32mcol_double()[39m,
psd_gam = [32mcol_double()[39m,
vol = [32mcol_double()[39m,
sizeclass = [31mcol_character()[39m,
lb = [32mcol_double()[39m,
ub = [32mcol_double()[39m,
binsize = [32mcol_double()[39m,
TotalParticles = [32mcol_double()[39m,
nparticles = [32mcol_double()[39m,
n_nparticles = [32mcol_double()[39m,
biovolume = [32mcol_double()[39m,
speed = [32mcol_double()[39m,
flux = [32mcol_double()[39m,
flux_fit = [32mcol_double()[39m,
GamPredictTP = [32mcol_double()[39m
)
uds <- read_csv("dataOut/unbinned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
PlotNParticles <- uds %>%
ggplot(aes(x = tot_nparticles, y = depth, col = profile)) +
facet_wrap(~project) +
geom_point(alpha = 0.3, shape = 1) +
scale_y_reverse() + scale_x_log10()
PlotNParticles
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
#FSG4 <- gam(tot_nparticles~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.002 0.106 84.96 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.204 1.366 3.346 0.0874 .
s(Day) 1.512 1.761 3.001 0.0386 *
s(Hour) 1.579 2.000 2.713 0.0389 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.183 Deviance explained = 24.2%
GCV = 0.73885 Scale est. = 0.67365 n = 60
#summary(FSG2)
#summary(FSG3)
#summary(FSG4)
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.07590024
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.05267319
summary(FSG3)$r.sq
[1] 0.0542546
But there is between projects:
ProjGam <- gam(tot_nparticles~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))
summary(ProjGam)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + factor(project)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.9983 0.4445 20.24 <2e-16 ***
factor(project)P16 17.0001 1.3754 12.36 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1 1 3.552 0.064 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.7 Deviance explained = 70.9%
GCV = 12.411 Scale est. = 11.855 n = 67
PlotPSDmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black"))
PlotParticlesmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
scale_x_log10() + theme(legend.position = "none")
PlotFluxmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = tot_flux_fit, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
scale_x_log10() + theme(legend.position = "none")
plot_grid(
PlotParticlesmany,
PlotPSDmany,
rel_widths = c(2, 3)
)
ggsave("figures/ParticlesPSDMany.png")
Saving 12 x 7.41 in image
ggsave("figures/ParticlesPSDMany.svg")
Saving 12 x 7.41 in image
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG4 <- gam(psd~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96083 0.01988 -199.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.713 1.917 56.197 5.5e-15 ***
s(Day) 1.000 1.000 1.078 0.3036
s(Hour) 1.604 2.000 3.049 0.0285 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.635 Deviance explained = 66.1%
GCV = 0.026006 Scale est. = 0.023702 n = 60
#summary(FSG2)
#summary(FSG3)
summary(FSG4)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96083 0.01995 -198.5 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.715 1.919 55.141 7.66e-15 ***
s(Hour) 1.565 2.000 2.497 0.0488 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.632 Deviance explained = 65.2%
GCV = 0.02572 Scale est. = 0.023885 n = 60
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.03804023
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] -0.004905911
summary(FSG3)$r.sq
[1] 0.6015425
Not a significant difference in PSD with respect to time.
But there is between projects:
ProjGam <- gam(psd~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))
summary(ProjGam)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + factor(project)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96167 0.02491 -159.029 <2e-16 ***
factor(project)P16 -0.19977 0.07708 -2.592 0.0118 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.224 1.398 34.54 2.26e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.437 Deviance explained = 45.6%
GCV = 0.039117 Scale est. = 0.037234 n = 67
I wonder if I can show that the profiles aren’t statistically significanlty different. Or that they are for that matter… I think in that case, I run a gam with and without a parameter for profile… And then quantify the effect size of that parameter
Or follow this Gavin Simpson Post https://fromthebottomoftheheap.net/2017/10/10/difference-splines-i/
or anova.gam {mgcv}
Calculate gams for each profile, and then run anova.gam to see if they are different…
PlotNParticlesEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
#geom_path(aes(x = tot_nparticles)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(2500, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
labs(x = "Particles/L", y = "depth (m)") +
theme(legend.position = "none") +
scale_shape_manual(values = c(1:5))
PlotNParticlesEP
I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)
PlotPSDEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
geom_path(aes(x = psd_gam)) +
geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(2500, 0)) + scale_color_manual(values = c("gray20", "brown")) +
scale_shape_manual(values = c(1:5)) + labs(y = "")
PlotPSDEP
I may just cow these togther.
plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))
Removed 611 rows containing missing values (geom_point).Removed 611 rows containing missing values (geom_point).Removed 611 row(s) containing missing values (geom_path).
ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")
Saving 10 x 4 in image
ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")
Saving 10 x 4 in image
mainParticleComponents <- bds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
select(project, profile, depth,
tot_nparticles, small_nparticles, big_nparticles,
tot_psd = psd, small_psd, big_psd,
tot_flux_fit, small_flux_fit, big_flux_fit) %>%
pivot_longer(cols = -c("project", "profile", "depth")) %>%
separate(name, c("size", "meas")) %>%
mutate(meas = recode(meas, nparticles = "particles/L")) %>%
mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))
Expected 2 pieces. Additional pieces discarded in 273 rows [7, 8, 9, 16, 17, 18, 25, 26, 27, 34, 35, 36, 43, 44, 45, 52, 53, 54, 61, 62, ...].
PlotFlx <- mainParticleComponents %>%
filter(meas != "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size ~ meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + theme(axis.title.x = element_blank(), legend.position = "none", strip.background.y = element_blank(), strip.text.y = element_blank(), plot.margin = unit(c(7,0,7,7), "pt")) + scale_color_manual(values = c("brown", "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = 175, color = "darkgreen")
PlotPSD <- mainParticleComponents %>%
filter(meas == "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size~meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = unit(c(7,7,26.5,0), "pt")) +
scale_color_manual(values = c("brown", "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = 175, color = "darkgreen")
plot_grid(PlotFlx, PlotPSD, rel_widths = c(3, 2))
Removed 246 rows containing missing values (geom_point).Removed 123 rows containing missing values (geom_point).
ggsave("figures/BigVsSmall.svg")
Saving 7.29 x 4.5 in image
ggsave("figures/BigVsSmall.png")
Saving 7.29 x 4.5 in image
Flux small and flux tot track so closely because ag > psd. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top.
eg_dataline <- bds %>%
filter(profile == "stn_043", depth == 162.5)
eg_slope = eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)
eg_datablock <- bes %>%
filter(profile == "stn_043", depth == 162.5)
eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)
eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)
EgNNP <- eg_datablock %>%
ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)")
EgNP <- eg_datablock %>%
ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Normalized Particles")
EgTP <- eg_datablock %>%
ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)")
plot_grid(EgNNP, EgTP, labels = c("A", "B"))
Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis
ggsave("figures/ExamplePSD163m.png")
Saving 7.29 x 4.5 in image
ggsave("figures/ExamplePSD163m.svg")
Saving 7.29 x 4.5 in image
bds %>%
ggplot(aes(y = depth, x = Flux_Smooth, col = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10()
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(Flux_Smooth~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(Flux_Smooth ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9334 37.21 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.875 1.984 11.823 0.000126 ***
s(Day) 1.756 1.939 2.168 0.094012 .
s(Hour) 1.114 2.000 1.658 0.059271 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.33 Deviance explained = 38.4%
GCV = 57.811 Scale est. = 52.277 n = 60
summary(FSG2)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9607 36.15 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.869 1.983 11.102 0.000206 ***
s(Day) 1.811 1.964 2.743 0.094614 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.29 Deviance explained = 33.4%
GCV = 60.062 Scale est. = 55.377 n = 60
summary(FSG3)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9944 34.93 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.857 1.979 10.73 0.000277 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.239 Deviance explained = 26.3%
GCV = 62.291 Scale est. = 59.325 n = 60
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.03975922
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.05063133
summary(FSG3)$r.sq
[1] 0.2392424
bds %>% filter(project == "ETNP") %>% select(profile, depth, Flux_Smooth) %>% pivot_wider(names_from = profile, values_from = Flux_Smooth)
Something is off. All of the flux profiles are identical. Skip this
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
plt1 <- bds %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = DFP, col = factor(time), shape = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + xlim(c(0.5, 1.5))+ geom_vline(xintercept = 1) +
scale_color_manual(values = c(rep("black", 5), rep("blue", 5))) + scale_shape_manual(values = rep(1:5, 2))
plotly::ggplotly(plt1)
What the heck is going on with DFP here. Why is it usually > 1 shouldn’t it be less than 1 when flux is decreasing? This very deep increasing flux seems improbable to me. Lets check the smooths. Or only go to 1000m.
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + facet_wrap(~project) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(2500, 0))+
scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d))) +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE)
pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlx
#plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse()+
#scale_x_log10() +
scale_x_log10(breaks = c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20)), limits = c(20, 180)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux") + theme(legend.position = "none")
pltFlxZoom
#plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time))) + facet_wrap(~project) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(2500, 0))+
scale_x_continuous(limits = c(-4.5, 1)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5})) + theme(legend.pos = "none")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3
#plotly::ggplotly(plt1pos)
#plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
Removed 15 rows containing missing values (geom_point).Removed 15 row(s) containing missing values (geom_path).
pgTop <- plot_grid(pltFlxNoLegend, pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
Removed 15 rows containing missing values (geom_point).Removed 15 row(s) containing missing values (geom_path).
pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""))
Removed 29 rows containing missing values (geom_point).Removed 29 row(s) containing missing values (geom_path).
pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
pgBoth
ggsave("figures/FluxDeepDive.png")
Saving 8 x 8 in image
ggsave("figures/FluxDeepDive.svg")
Saving 8 x 8 in image
#plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
Removed 15 rows containing missing values (geom_point).Removed 15 row(s) containing missing values (geom_path).
pgTop <- plot_grid(pltFlxNoLegend + ylim(c(1000, 0)), pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Removed 42 rows containing missing values (geom_point).Removed 42 row(s) containing missing values (geom_path).
pgBottom <- plot_grid(pltDelta3 + ylim(c(1000, 0)), pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""))
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Removed 56 rows containing missing values (geom_point).Removed 56 row(s) containing missing values (geom_path).
pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
pgBoth
ggsave("figures/FluxShallowDive.png")
Saving 8 x 8 in image
ggsave("figures/FluxShallowDive.svg")
Saving 8 x 8 in image
Test for day to day and hourly variability
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
DFG1 <- gam(DF/DZ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
DFG2 <- gam(DF/DZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
DFG3 <- gam(DF/DZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(DFG1)
Family: gaussian
Link function: identity
Formula:
DF/DZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1001 0.0171 -5.854 2.94e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.911 1.992 22.994 7.83e-08 ***
s(Day) 1.870 1.983 7.097 0.00126 **
s(Hour) 1.292 2.000 2.013 0.05235 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.514 Deviance explained = 55.5%
GCV = 0.019516 Scale est. = 0.017541 n = 60
summary(DFG2)
Family: gaussian
Link function: identity
Formula:
DF/DZ ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.10008 0.01769 -5.657 5.67e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.902 1.990 22.110 1.37e-07 ***
s(Day) 1.904 1.991 6.768 0.00339 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.479 Deviance explained = 51.3%
GCV = 0.020417 Scale est. = 0.018782 n = 60
summary(DFG3)
Family: gaussian
Link function: identity
Formula:
DF/DZ ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.10008 0.01944 -5.149 3.37e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.878 1.985 19.24 9.12e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.371 Deviance explained = 39.1%
GCV = 0.023812 Scale est. = 0.02267 n = 60
summary(DFG1)$r.sq - summary(DFG2)$r.sq
[1] 0.03441623
summary(DFG2)$r.sq - summary(DFG3)$r.sq
[1] 0.1078441
summary(DFG3)$r.sq
[1] 0.3713016
plot(DFG1)
#Osps
(u mol C / m^3 / day)
bds %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + facet_wrap(~project) + geom_point(size = 2) + scale_y_reverse(limits = c(2500, 0)) + geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = bquote("Flux Size Shift"~μmol/m^3/day)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue")
#ggsave("..figures/FluxSizeShift.svg"
ggsave("figures/FluxSizeShift.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FluxSizeShift.svg")
Saving 7.29 x 4.5 in image
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(OZG1)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.024817 0.002245 11.05 1.48e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.8892 1.988 8.422 0.000434 ***
s(Day) 1.8718 1.983 6.012 0.003119 **
s(Hour) 0.5513 2.000 0.396 0.243644
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.347 Deviance explained = 39.5%
GCV = 0.00033186 Scale est. = 0.00030248 n = 60
summary(OZG2)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.02482 0.00226 10.98 1.67e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.89 1.988 8.429 0.000428 ***
s(Day) 1.89 1.988 6.748 0.003790 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.338 Deviance explained = 38.1%
GCV = 0.00033306 Scale est. = 0.00030653 n = 60
summary(OZG3)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.024817 0.002479 10.01 3.55e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.872 1.984 7.338 0.000983 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.204 Deviance explained = 22.9%
GCV = 0.00038738 Scale est. = 0.00036884 n = 60
summary(OZG1)$r.sq - summary(OZG2)$r.sq
[1] 0.008753478
summary(OZG2)$r.sq - summary(OZG3)$r.sq
[1] 0.1344912
summary(OZG3)$r.sq
[1] 0.2038316
bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = ospsDZ)) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(500, 0)) + geom_vline(xintercept = 0)
trapFlux3 <- read_csv("dataOut/fluxMS_distilled.csv")
Parsed with column specification:
cols(
Class = [31mcol_character()[39m,
Depth = [32mcol_double()[39m,
TrapID = [31mcol_character()[39m,
TrapType = [31mcol_character()[39m,
SampleType = [31mcol_character()[39m,
C_flux = [32mcol_double()[39m,
C_flux_umol = [32mcol_double()[39m
)
UVPFluxComb <- read_csv("dataOut/CombinedProfileFluxEst_DS.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
Flux = [32mcol_double()[39m
)
trapFlux3
UVPFluxComb
fluxMS_distilled_toPlot <- trapFlux3 %>%
mutate(SampleType = recode(SampleType, `plus.p` = "plus-particles", top = "top-collector"))
UVPFluxComb %>%
ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(0, 200)) +
geom_point(aes(y = Depth, x = C_flux_umol, fill = SampleType, shape = TrapType),
colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot) +
geom_point(aes(x = Flux), size = 3, shape = 21, color = "white", fill = "black") +
geom_point(aes(x = -1, y = -1, size = "UVP")) + # dummy point for the legend
scale_shape_manual(values = c(25, 22))+
scale_size_manual(values = 1, name = "") +
ylab("Depth (m)") + xlab("Flux µmolC/m^2/day") +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
theme_cowplot() +
theme(
legend.position = c(0.5, 0.4),
legend.box.background = element_rect(color = "black", size = 0.5),
legend.margin = margin(-10, 5, 10, 5)
) +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
ggsave("figures/FittedFlux.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FittedFlux.svg")
Saving 7.29 x 4.5 in image
SameGam <- gam(TotalParticles ~s(log(lb), log(depth), by = factor(time)), offset = log(vol * binsize), family = nb(),
data = bes %>% filter(project == "ETNP"))
Pred <- exp(predict(SameGam, Expanded))
Pred <- exp(predict(SameGam, Expanded))
meanBese <- ToPlot %>% filter(lb <= 2) %>% group_by(lb, depth) %>% summarize(nparticles = mean(nparticles), nnparticles = mean(nnparticles))
WBColorMap <- meanBese%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c(name = "log10(number density \n (normalized))") + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue")
WBColorMap
Average of everything
meanBese043 <- ToPlot %>% filter(lb <= 2, time == "2017-01-13 11:51:31")
meanBese043%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c() + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 175, color = "darkgreen")
Just 043
mbGam <- meanBese %>% group_by(depth) %>% nest() %>%
mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>%
mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue")
mbGam <- meanBese043 %>% group_by(depth) %>% nest() %>%
mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>%
mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
pWBPSD <- mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue")
pWBPSD
bds %>% filter(profile == “stn_043”, depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = “darkgreen”) + geom_hline(yintercept = 950, color = “darkblue”) + geom_ribbon(alpha = 0.2) + labs(x = “PSD slope”)
All of them
bds %>% filter(profile == "stn_043", depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue") + geom_ribbon(alpha = 0.2) + labs(x = "PSD slope")
043 only
bds %>% filter(profile == "stn_043", depth <= 2000, depth > 175) %>% ggplot(aes(x = small_biovolume, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue") + geom_point()
ubDf0 <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global)
ubDf <- ubDf0 %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth)
photicBiomass <- ubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
ubDf <- ubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length
ubDf %>% ggplot(aes(x = nbiomass, y = depth , group = time, col = time)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1))
ubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- ubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
ubDf <- ubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
ubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 175, color = "darkgreen")
PubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb < 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- PubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
PubDf <- PubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBS <- PubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1.2)) + geom_hline(yintercept = 175, color = "darkgreen") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = 950, color = "darkblue") + labs( x = "Small particle mass (norm.)")
pWBS
LubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb >= 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- LubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
LubDf <- LubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBL <- LubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 175, color = "darkgreen") + labs( x = "Large particle mass (norm.)") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = 950, color = "darkblue")
pWBL
For tom and danielle
WBColorMap
pWBPSD
pWBS
pWBL
WBFig5 <- plot_grid(pWBPSD, pWBS,pWBL, nrow = 1, labels = c("B", "C", "D"))
Removed 6 row(s) containing missing values (geom_path).Removed 8 row(s) containing missing values (geom_path).
WBFig5
WBcombined <- plot_grid(WBColorMap + theme(plot.margin = unit(c(0,3,0, 3), "cm")), WBFig5, ncol = 1, labels = c("A", ""))
WBcombined
ggsave("figures/WBModelValidation.png")
Saving 8 x 6 in image
dataBinned <- read_csv("data/backscatter_table_go7.csv")
dataBinned_01 <- dataBinned %>%
mutate(timeMex = with_tz(time_bin, tzone = "US/Central") )
startDay <- dataBinned_01$timeMex %>% na.omit %>% min %>% floor_date(unit = "days")
endDay <- dataBinned_01$timeMex %>% na.omit %>% max %>% ceiling_date(unit = "days")
timeBreaks <- seq(from = startDay, to = endDay, by = "12 hours")
timeLabels <- format(timeBreaks)
dataBinned_01 %>% filter(frequency == 18000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -75), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
#ggsave("figures/stationP2_EK60_18kOnly.svg")
ggsave("figures/stationP2_EK60_18kOnly.png")
dataBinned_01 %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -75), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
facet_wrap(~frequency, ncol = 1)
#ggsave("figures/stationP2_EK60_go7.svg", width = 8, height = 20)
ggsave("figures/stationP2_EK60_go7.png", width = 8, height = 20)